МІНІСТЕРСТВО ОСВІТИ І НАУКИ УКРАЇНИ
НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ “ЛЬВІВСЬКА ПОЛІТЕХНІКА”
РОЗРАХУНКОВА РОБОТА №2
З дисципліни : Алгоритмізація і програмування.
Знаходження динамічних характеристик технологічного об’єкта шляхом розв’язування системи диференціальних рівнянь.
Варіант №2
Завдання
Для протічної гідравлічної ємності з коротким та довгим трубопроводами відповідно на вході та виході математична модель має вигляд :
де
Знайти розв’язок системи звичайних диференціальних рівнянь (ЗДР) – H(t), Q1(t) на інтервалі t є [0; 150] с. з кроком Δt=0.5 c, якщо
d=1 м
ξ=0.9
ρ=1000 кг/м³
P1=50000 Па
L1=95 м
r1=0.12 м
Q=0.02 м³/c
ν=0.00001 Па*с
Початкові умови : значення рівня H(0)=1 м та витрати Q(0)=0,02 м³/c.
Розв’язання
Розв’яжемо систему відносно перших похідних функцій, тобто
Позначимо
,
Тоді схема числового інтегрування системи двох ЗДР матиме вигляд :
Блок-схема методу Рунге-Кута для розв’язання системи двох ЗДР
Програма чисельного роз’язування системи двох ЗДР
#include<stdio.h>
#include<conio.h>
#include<math.h>
float f1(float t, float H, float Q1)
{float dz=0.9,po=1000,v=1e-5,P1=50000,L1=95,r1=0.12,d=1,Q=0.02,g=9.81;
float A=4*M_PI*pow(r1,3)/dz;
float k1=2*M_PI*r1*r1*sqrt(r1/(L1*dz));
return (1/A)*(k1*k1*(P1-po*g*H)/po-Q1*Q1);
}
float f2(float t, float H, float Q1)
{float Q=0.02, d=1;
float S=M_PI*d*d/4;
return (Q1-Q)/S;
}
main()
{clrscr();
float t,t0=0,tf=150,n=300;
float h,k1,k2,k3,k4,l1,l2,l3,l4,Hi=1,Q1i=0.02;
h=(tf-t0)/n;
for(t=t0;t<=tf;t+=h)
{printf(«\n t=%f Hi=%f Q1i=%f»,t,Hi,Q1i);
k1=h*f1(t,Hi,Q1i);
l1=h*f2(t,Hi,Q1i);
k2=h*f1(t+h/2,Hi+k1/2,Q1i+l1/2);
l2=h*f2(t+h/2,Hi+k1/2,Q1i+l1/2);
k3=h*f1(t+h/2,Hi+k2/2,Q1i+l2/2);
l3=h*f2(t+h/2,Hi+k2/2,Q1i+l2/2);
k4=h*f1(t+h,Hi+k3,Q1i+l3);
l4=h*f2(t+h,Hi+k3,Q1i+l3);
Hi=Hi+(k1+2*k2+2*k3+k4)/6;
Q1i=Q1i+(l1+2*l2+2*l3+l4)/6;
}
getch();
return 0;
}
Запишемо в таблицю кожне тридцяте значення функцій для рівня H(t), м та витрати Q1(t), м³/с.
t,c
0
15
30
45
60
75
H(t),м
1.0000
1.0370
1.0716
1.1038
1.1339
1.1619
Q1(t),м³/с.
0.0200
0.0200
0.0200
0.0200
0.0200
0.0200
t,c
90
105
120
135
150
H(t),м
1.1880
1.2124
1.2351
1.2563
1.2760
Q1(t),м³/с.
0.0200
0.0200
0.0200
0.0200
0.0200
Розв’язок системи двох ЗДР засобами MATLAB
Створимо файл data.m в якому запишемо вхідні дані та обчислимо коефіцієнти
%data
dz=0.9;po=1000;v=1e-5;g=9.81;P1=50000;
L1=95;r1=0.12;d=1;Q=0.02;
A=(4*pi*r1^3)/dz;
K1=2*pi*r1^2*sqrt(r1/(L1*dz));
S=pi*d^2/4;
Файл-функція ff.m, в якому записані праві частини диференціальних рівнянь, розв’язаних відносно перших похідних, має вигляд :
function y=ff(t,x);
data;
H=x(1);Q1=x(2);
y=[(1/A)*((K1^2/po)*(P1-po*g*H)-Q1^2);
(1/S)*(Q1-Q)];
Дані та обчислені значення коефіцієнтів передають в файл-функцію ff.m шляхом включення файлу даних data.m.
Файл main.m, в якому задають інтервал часу (t0, tf), на якому інтегрують систему ЗДР, початкові умови (x0=[1 0.02]) та записують функцію MATLAB ode45(), призначену для розв’язку ЗДР чи системи ЗДР числовим методом має вигляд :
%main
t0=0;tf=150;x0=[1 0.02];
[t,y]=ode45('ff',[t0 tf],x0)
figure(1);subplot(2,1,1);plot(t,y(:,1),'k');grid;xlabel('t,c');ylabel('H,m');
subplot(2,1,2);plot(t,y(:,2),'k');grid;xlabel('t,c');ylabel('Q_1,m^3/c');
Графіки перехідних прцесів – зміни рівня в ємності та витрати рідини в трубопроводі.
Накладемо графіки, отримані в MATLAB та значення отримані в С. Для цього модефікуємо файл main.m, додавши до нього дані з таблиці :
%main
t0=0;tf=150;x0=[1 0.02];
[t,y]=ode45('ff',[t0 tf],x0)
ti=[0:15:150];
Hi=[1 1.0370 1.0716 1.1038 1.1339 1.1619 1.1880 1.2124 1.2351 1.2563 1.2760];
Q1i=[0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02];
figure(1);subplot(2,1,1);plot(t,y(:,1),'k',ti,Hi,'b*');grid;xlabel('t,c');ylabel('H,m');
subplot...
